# Vaccine Diplomacy

## Replicates figures in main text and appendix

# Figure 1: Cumulative doses per 100 people across six Latin American countries and survey dates ----

df <- read.csv("raw vaccination data.csv")
df <- df |> filter(location %in% c('Argentina', 'Brazil', 'Chile', 'Colombia', 'Mexico', 'Peru'))
df$location <- factor(df$location)
df <- df |> select(location, date, total_vaccinations_per_hundred)

df$date <- as.character(df$date)
df$date <- as.Date(df$date, "%m/%d/%Y")

df <- df  |>  filter(date <= as.Date("2021-06-30"))

df <- df |> mutate(location = case_when(
  location=="Argentina" ~ "ARG", 
  location=="Brazil" ~ "BRA", 
  location=="Chile" ~ "CHI", 
  location=="Colombia" ~ "COL", 
  location=="Mexico" ~ "MEX", 
  location=="Peru" ~ "PER"
))

rollout <- df |> 
  ggplot(aes(x=date, y=total_vaccinations_per_hundred, group=location, color=location)) +
  geom_textline(aes(linetype = location, label = location),
                data = df[!is.na(df$total_vaccinations_per_hundred), ],
                linewidth = 1,
                size = 2.75,
                hjust = 0.875) +
  theme_classic()+
  scale_fill_grey() + 
  xlab("") + 
  ylab("Cumulative doses administered per 100 people")  +
  scale_color_manual(values=c('#000000','#000000', '#000000','#949494', '#949494', '#949494')) +
  scale_linetype_manual(values=c("solid", "twodash", "longdash", "solid", "dashed", "dotdash")) +
  theme(legend.title = element_blank(),
        legend.text = element_blank(),
        legend.position = "none",
        axis.text.x = element_text(size=15),
        axis.text.y = element_text(size=15),
        axis.title.y = element_text(size=15)) +
  scale_x_date(date_labels="%b %y",date_breaks  ="1 month") +
  geom_label(label = "Baseline",
             x = as.numeric(as.Date("2021-01-21")),
             y = 120, label.size = NA, color='black', 
             size=3.5) +
  geom_label(label = "Endline",
             x = as.numeric(as.Date("2021-06-07")),
             y = 120, label.size = NA, color='black', 
             size=3.5) +
  geom_vline(xintercept = as.numeric(as.Date("2021-01-11")), linetype=3, size=1) +
  geom_vline(xintercept = as.numeric(as.Date("2021-01-30")), linetype=3, size=1) +
  geom_vline(xintercept = as.numeric(as.Date("2021-05-21")), linetype=3, size=1) +
  geom_vline(xintercept = as.numeric(as.Date("2021-06-24")), linetype=3, size=1)

plot(rollout)

ggsave(filename = "Figures/rollout.pdf",
       width = 1.375*5.5,
       height = 1*5.5,
       units = "in",
       dpi = 320)

rm(df, rollout)

# Figure 2: Overview of endline survey flow and treatment assignments ----

## Figure made in LaTeX

# Figure 3: Number of vaccine doses per adult from each vaccine-developer country (in May 2021) ----

doses_received <- read.csv("doses_received.csv")

doses_received <- doses_received |> 
  mutate(adult_pop = case_when(
    receiving_country == "Argentina" ~ 31541111,
    receiving_country == "Brasil" ~ 167709191,
    receiving_country == "Chile" ~ 12757617,
    receiving_country == "Colombia" ~ 39563361,
    receiving_country == "Mexico" ~ 82997882,
    receiving_country == "Peru" ~ 23064256))

diplomacy_wide <- read.csv("vaccine_diplomacy_wide.csv",
                           encoding = "UTF-8")

diplomacy_wide <- diplomacy_wide |> 
  mutate(know_correct_vax = case_when(
    vax_which_end %in% c("CanSino", "Sinopharm", "Sinovac (CoronaVac)") & vax_know_orig_end == "China" ~ 1,
    vax_which_end %in% c("Covishield", "Covaxin") & vax_know_orig_end %in% c("India", "?\u008dndia") ~ 1,
    vax_which_end %in% c("Sputnik V") & vax_know_orig_end %in% c("Rusia", "Russia") ~ 1,
    vax_which_end %in% c("Astra-Zeneca/Oxford") & vax_know_orig_end %in% c("Reino Unido", "Reino Unido/Gra Bretanha") ~ 1, 
    vax_which_end %in% c("Johnson and Johnson (Janssen)", "Moderna", "Pfizer/BioNTech") & vax_know_orig_end == "Estados Unidos" ~ 1, 
    ),
    vax_country = factor(case_when(
      vax_which_end %in% c("CanSino", "Sinopharm", "Sinovac (CoronaVac)") ~ "China",
      vax_which_end %in% c("Covishield", "Covaxin") ~ "India",
      vax_which_end %in% c("Sputnik V") ~ "Russia",
      vax_which_end %in% c("Astra-Zeneca/Oxford") ~ "United Kingdom",
      vax_which_end %in% c("Johnson and Johnson (Janssen)", "Moderna", "Pfizer/BioNTech") ~ "United States"
    )),
    doses = case_when(
      vax_doses_end %in% c("Dois (2)", "Dos (2)") ~ 2,
      vax_doses_end %in% c("Um (1)", "Una (1)") ~ 1,
      TRUE ~ 0
    ))

diplomacy <- diplomacy_wide

rm(diplomacy_wide)

country_data <- lapply(sort(unique(diplomacy$country)), function(i)
  filter(diplomacy, country == i)  
)

country_plots_survey <- lapply(country_data, function(i)
 aggregate(i$doses, by = list(i$vax_country), FUN = "sum") |> 
   rename(vax_country = Group.1,
          doses = x) |> 
   mutate(total_ppl = length(unique(i$mergeid)),
          source = "Survey Pool")
)

country_plots_nd <- lapply(sort(unique(diplomacy$country)), function(i)
  filter(doses_received, receiving_country == i) |> 
    select(-c(receiving_country, pct_total)) |> 
    mutate(source = "National Data") |> 
    rename(vax_country = sending_country,
           total_ppl = adult_pop)
)

country_plots <- lapply(1:length(unique(diplomacy$country)), function(i)
   rbind(country_plots_survey[[i]], country_plots_nd[[i]])
)

rm(country_plots_survey, country_plots_nd)

country_plot <- lapply(1:length(unique(diplomacy$country)), function(i)
  data.frame(vax_country = rep(unique(diplomacy$vax_country), 2),
             source = c(rep("Survey Pool", length(unique(diplomacy$vax_country))),
                        rep("National Data", length(unique(diplomacy$vax_country))))) |> 
    na.omit()
)

for (i in 1:length(country_plots)) {
  
  country_plot[[i]] <- left_join(country_plot[[i]], country_plots[[i]]) |> 
    mutate(doses = ifelse(is.na(doses), 0, doses),
           per_cap = doses/total_ppl,
           total_ppl = ifelse(is.na(total_ppl), 0, total_ppl),
           per_cap = ifelse(is.na(per_cap), 0, per_cap),
           vax_country = case_when(
             vax_country == "United Kingdom" ~ "UK", 
             vax_country == "United States" ~ "US", 
             vax_country == "Russia" ~ "Russia", 
             vax_country == "China" ~ "China", 
             vax_country == "India" ~ "India", 
             vax_country == "UK" ~ "UK",
             vax_country == "US" ~ "US" 
           ))
  
}

rm(country_plots, i)

names(country_plot) <- c("Argentina", "Brazil", "Chile", "Colombia", "Mexico", "Peru")

country_barplot <- lapply(names(country_plot), function(i)
  
  ggplot(data = country_plot[[i]], aes(x = vax_country, y = per_cap, fill = source)) +
    geom_bar(position="dodge", color = 'black', stat="identity") +
    ylab("") +
    xlab("") +
    ggtitle(i) +
    theme(plot.title = element_text(size = 23, hjust = 0.5)) + 
    theme(legend.position="bottom", legend.title = element_blank()) +
    coord_cartesian(ylim=c(0, 1.25)) + 
    scale_fill_brewer(palette="Blues", labels = c("National Data", "Survey Pool")) + 
    theme(axis.title = element_text(size = 15), 
          legend.text = element_text(size = 20), 
          axis.text = element_text(size=20),
          legend.margin = margin(t = -15))
)

names(country_barplot) <- paste0(tolower(names(country_plot)), "_doses_barplot.pdf")

country_barplot

for (i in names(country_barplot)) {
  
  ggsave(filename = paste0("Figures/", i),
         plot = country_barplot[[i]],
         width = 1*5,
         height = 1*5,
         units = "in",
         dpi = 320)
  
}

rm(country_barplot,
   country_data,
   country_plot,
   diplomacy,
   doses_received,
   i)

# Figure 4: Respondents perceptions of motivation for developing vaccines, by vaccine developing country ----

diplomacy_wide <- read.csv("vaccine_diplomacy_wide.csv",
                           encoding = "UTF-8")

mech_df <- diplomacy_wide |> 
  select(mech_china_stopcv, mech_china_help, mech_china_supp, mech_china_dep, mech_china_econ,
         mech_india_stopcv, mech_india_help, mech_india_supp, mech_india_dep, mech_india_econ,
         mech_russia_stopcv, mech_russia_help, mech_russia_supp, mech_russia_dep, mech_russia_econ,
         mech_uk_stopcv, mech_uk_help, mech_uk_supp, mech_uk_dep, mech_uk_econ,
         mech_us_stopcv, mech_us_help, mech_us_supp, mech_us_dep, mech_us_econ)

mech_df <- mech_df |> pivot_longer(mech_china_stopcv:mech_us_econ,
                     names_to = "mechanism",
                     values_to = "count") |>
  filter(count == 1) |> 
  mutate(country = str_replace(mechanism, "mech_", ""),
         country = str_replace(country, "_.*$", ""),
         country = case_when(
           country == "china" ~ "China",
           country == "india" ~ "India",
           country == "russia" ~ "Russia",
           country == "uk" ~ "United Kingdom",
           country == "us" ~ "United States"
         ),
         mechanism = str_replace(mechanism, "mech_.*_", ""),
         mechanism = case_when(
             mechanism == "stopcv" ~ "Stop Covid-19 Spread", 
             mechanism == "help" ~ "Help Country", 
             mechanism == "supp" ~ "Increase Support", 
             mechanism == "dep" ~ "Increase Dependence", 
             mechanism == "econ" ~ "Obtain Profits"
         ),
         mechanism = factor(mechanism, levels=c("Stop Covid-19 Spread",
                                       "Help Country", 
                                       "Increase Support", 
                                       "Increase Dependence", 
                                       "Obtain Profits"))) |> 
  group_by(country, mechanism)|> 
  summarise(mechanism_count = sum(count)) |> 
  mutate(mechanism_prop = mechanism_count/sum(mechanism_count))

mech_barplot <- lapply(sort(unique(mech_df$country)), function(i)
  
  mech_df |> 
    filter(country == i) |> 
  ggplot(aes(x=factor(mechanism), y = mechanism_prop, fill=factor(mechanism))) +
    geom_bar(stat = "identity", color = 'black', fill = 'lightblue') +
    ylab("") + 
    xlab("") +
    ggtitle(i) +
    scale_y_continuous(labels = scales::percent_format(accuracy = 5L)) +
    scale_x_discrete(labels = c("Stop\nCOVID-19\nspread",
                                "Help\ncountry", 
                                "Increase\nsupport", 
                                "Increase\ndependence", 
                                "Obtain\nprofits")) +
    theme(plot.title = element_text(size = 20, hjust = 0.5)) +
    theme(axis.text = element_text(size=12)) +
    coord_cartesian(ylim=c(0, .5)) 
  
)

names(mech_barplot) <- paste0(c("china",
                                "india",
                                "russia",
                                "uk",
                                "us"), "_mech.pdf")

mech_barplot

for (i in names(mech_barplot)) {
  
  ggsave(filename = paste0("Figures/", i),
         plot = mech_barplot[[i]],
         width = 1*5,
         height = 1*5,
         units = "in",
         dpi = 320)
  
}

rm(diplomacy_wide, mech_barplot, mech_df, i)

# Figure 5: Perceived and actual ranking of doses delivered by each vaccine developer country in each country surveyed ----

heatmap_df <- read.csv("heatmap_data.csv")

heatmap_df <- heatmap_df |> 
  mutate(Sending.Country = str_replace(Sending.Country, " \\(Tied", "\\\n(Tied"))

arg <- heatmap_df |> filter(Country=="Argentina")
arg$Sending.Country <- factor(arg$Sending.Country, levels = c("US\n(Tied 4th)",
                                                              "UK\n(Tied 4th)",
                                                              "India (3rd)",
                                                              "China (2nd)",
                                                              "Russia (1st)"))
arg$Perceived.Rank <- factor(arg$Perceived.Rank, levels = c("5", "4", "3", "2", "1"))

br <- heatmap_df |> filter(Country=="Brazil")
br$Sending.Country <- factor(br$Sending.Country, levels = c("Russia\n(Tied 4th)", "India\n(Tied 4th)",  "US (3rd)",  "UK (2nd)", "China (1st)"))
br$Perceived.Rank <- factor(br$Perceived.Rank, levels = c("5", "4", "3", "2", "1"))

cl <- heatmap_df |> filter(Country=="Chile")
cl$Sending.Country <- factor(cl$Sending.Country, levels = c("India\n(Tied 4th)", "Russia\n(Tied 4th)", "UK (3rd)",  "US (2nd)", "China (1st)"))
cl$Perceived.Rank <- factor(cl$Perceived.Rank, levels = c("5", "4", "3", "2", "1"))

co <- heatmap_df |> filter(Country=="Colombia")
co$Sending.Country <- factor(co$Sending.Country, levels = c("India\n(Tied 4th)", "Russia\n(Tied 4th)","UK (3rd)",  "US (2nd)", "China (1st)"))
co$Perceived.Rank <- factor(co$Perceived.Rank, levels = c("5", "4", "3", "2", "1"))

mx <- heatmap_df |> filter(Country=="Mexico")
mx$Sending.Country <- factor(mx$Sending.Country, levels = c("India (5th)", "Russia (4th)", "UK (3rd)", "China (2nd)", "US (1st)"))
mx$Perceived.Rank <- factor(mx$Perceived.Rank, levels = c("5", "4", "3", "2", "1"))

pe <- heatmap_df |> filter(Country=="Peru")
pe$Sending.Country <- factor(pe$Sending.Country, levels = c("India\n(Tied 3rd)", "Russia\n(Tied 3rd)", "UK\n(Tied 3rd)", 
                                                            "China (2nd)", 
                                                            "US (1st)"))
pe$Perceived.Rank <- factor(pe$Perceived.Rank, levels = c("5", "4", "3", "2", "1"))

heatmap_list <- list("Argentina" = arg,
                     "Brazil" = br,
                     "Chile" = cl,
                     "Colombia" = co,
                     "Mexico" = mx,
                     "Peru" = pe)

rm(arg, br, cl, co, mx, pe)

heatmap_figs <- lapply(1:length(heatmap_list), function(i)
  
  heatmap_list[[i]] |> 
    ggplot(aes(as.factor(Perceived.Rank), as.factor(Sending.Country), fill=X)) +
    geom_tile(color="black") +
    geom_text(data=heatmap_list[[i]],
              aes(label=Percentage), size=4.2, fontface="bold") +
    scale_fill_gradient(low="white", high="#055795") +
    theme(panel.background = element_rect(fill = 'white'), 
          axis.ticks.x=element_blank(),
          axis.ticks.y=element_blank(), 
          plot.title = element_text(size = 18, face = "bold")
    ) +
    theme(legend.position = "none", 
          axis.text=element_text(size=12),
          axis.title=element_text(size=14)) +
    coord_fixed() +
    ylab("Actual Rank of Doses Delivered by Country") +
    xlab("Perceived Rank") +
    ggtitle(paste0("Respondents in ", names(heatmap_list)[[i]]))
  
)

names(heatmap_figs) <- paste0("heatmap_", c("arg",
                                            "br",
                                            "cl",
                                            "co",
                                            "mx",
                                            "pe"), ".pdf")

heatmap_figs

for (i in names(heatmap_figs)) {
  
  ggsave(filename = paste0("Figures/", i),
         plot = heatmap_figs[[i]],
         width = 1*5,
         height = 1*5,
         units = "in",
         dpi = 320)
  
}

rm(heatmap_df, heatmap_figs, heatmap_list, i)

# Figure 6: Screenshot of the information treatment (from Argentina) ----

## Screenshot from Qualtrics

# Figure 7: Moderation of the effect of the aggregate vaccine distribution information treatment on trust in foreign governments, by vaccine-developer country rank ----

diplomacy_long <-
  read.csv("vaccine_diplomacy_long.csv", encoding = "UTF-8")

# Marginal effects plot

marginal_effects_models <- function(outcome, interaction) {
  lm_robust(
    as.formula(paste0(outcome, "~ treat_rank*", interaction, " + factor(country_pre_trust)*factor(country_exp)")),
    fixed_effects = exp_block,
    clusters = mergeid,
    se_type = "stata",
    data = diplomacy_long
  )
}

rank_model <- marginal_effects_models(outcome = "country_post_trust",
                                      interaction = "factor(rev_rank_reported_ties)")

rank_model_binary <- marginal_effects_models(outcome = "country_post_trust_binary",
                                      interaction = "factor(rev_rank_reported_ties)")

## partial Y / partial T = b + d * M

marginal_effect <- function(model, treatment_name, moderator_name) {
  
  result <- c(model$coefficients[treatment_name],
    model$coefficients[treatment_name] +
      model$coefficients[grep(paste0(treatment_name, ":"), names(model$coefficients))])
  
  return(as.numeric(result))
  
}

marginal_effect_bound <- function(model, treatment_name, moderator_name, bound) {
  
  marginal_effect <- c(model$coefficients[treatment_name],
                       model$coefficients[treatment_name] +
                         model$coefficients[grep(paste0(treatment_name, ":"), names(model$coefficients))])
  
  ## V( partial Y / partial T ) = V(b) + M^2 * V(d) + 2 * M * Cov(b,d),
  
  vcov_matrix <- data.frame(model$vcov)
  
  M_indices <- grep(paste0(treatment_name, "\\."), names(vcov_matrix))
  
  v_partialy_partialt <-
    c(vcov_matrix[treatment_name, treatment_name],
      sapply(1:length(M_indices), function(i)
        vcov_matrix[treatment_name, treatment_name] +
          vcov_matrix[M_indices[i], M_indices[i]] +
          2 * vcov_matrix[treatment_name, M_indices[i]]
      )
    )
  
  v_partialy_partialt <- sqrt(v_partialy_partialt)
  
  if (bound == "upper") {
    result <- marginal_effect + (1.96*v_partialy_partialt)
  } else if (bound == "lower") {
    result <- marginal_effect - (1.96*v_partialy_partialt)
  }
  
  return(as.numeric(result))
  
}

# Marginal effects plots for rank

rank_plots_df <-
  lapply(list(rank_model, rank_model_binary), function(i)
    cbind.data.frame(
      m = sort(unique(diplomacy_long$rev_rank_reported_ties)),
      m_effect = marginal_effect(
        model = i,
        treatment_name = "treat_rank",
        moderator_name = "factor(rev_rank_reported_ties)"
      ),
      m_lower = marginal_effect_bound(
        model = i,
        treatment_name = "treat_rank",
        moderator_name = "factor(rev_rank_reported_ties)",
        bound = "lower"
      ),
      m_upper = marginal_effect_bound(
        model = i,
        treatment_name = "treat_rank",
        moderator_name = "factor(rev_rank_reported_ties)",
        bound = "upper"
      )
    ))

rank_plots <- lapply(rank_plots_df, function(i)
  ggplot(i,
       aes(x = m,
           y = m_effect)) +
  geom_point() +
  geom_errorbar(aes(x = m, ymin = m_lower, ymax = m_upper, width = 0)) +
  geom_hline(yintercept = 0, linetype = "dotted") +
  scale_y_continuous("Marginal effect of treatment") +
  theme(axis.line = element_line(size = 0.5),
        axis.title.y = element_text(size = 18),
        axis.text.y = element_text(size = 15),
        axis.ticks.x = element_blank(), 
        axis.text.x =element_blank(),
        axis.title.x = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        plot.margin = unit(c(0, 0, -0.2, 0), "cm"))
)

names(rank_plots) <- c("scale", "binary")

rank_plot_density <- ggplot(diplomacy_long, aes(x = rev_rank_reported_ties)) +
  geom_histogram(fill = "grey", bins = 9, binwidth = 0.25) +
  scale_x_continuous("Reversed rank") +
  theme(axis.line = element_line(size = 0.5),
        axis.title.x = element_text(size = 18),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(), 
        axis.title.y = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        plot.margin = unit(c(0, 0, 0, 0), "cm"))

ggsave(filename = "Figures/rank_marginaleffects.pdf",
       plot = plot_grid(rank_plots[["scale"]], rank_plot_density, align = "v", nrow = 2, rel_heights = c(5/6, 1/6)),
       width = 1*5,
       height = 1*5,
       units = "in",
       dpi = 320)

ggsave(filename = "Figures/rank_marginaleffects_binary.pdf",
       plot = plot_grid(rank_plots[["binary"]], rank_plot_density, align = "v", nrow = 2, rel_heights = c(5/6, 1/6)),
       width = 1*5,
       height = 1*5,
       units = "in",
       dpi = 320)

# Figure 8: Moderation of the effect of the aggregate vaccine distribution information on the perceived motivation of government of the country where the vaccine was developed for distributing vaccines, by vaccine-developer country rank ----

mechanism_outcomes <- c(
  "Stop COVID-19 spread"          = "country_post_mechanisms_stopcovid",
  "Help respondent country"       = "country_post_mechanisms_help",
  "Increase support for sender"   = "country_post_mechanisms_support",
  "Increase dependence on sender" = "country_post_mechanisms_dependence",
  "Obtain economic profits"       = "country_post_mechanisms_economics"
)

mechanism_models <- lapply(mechanism_outcomes, function(y)
  marginal_effects_models(outcome = y,
                          interaction = "factor(rev_rank_reported_ties)")
)

mechanism_plots_df <-
  lapply(mechanism_models, function(i)
    cbind.data.frame(
      m = as.character(c(2, 3, 4, 5)),
      m_effect = marginal_effect(
        model = i,
        treatment_name = "treat_rank",
        moderator_name = "factor(rev_rank_reported_ties)"
      ),
      m_lower = marginal_effect_bound(
        model = i,
        treatment_name = "treat_rank",
        moderator_name = "factor(rev_rank_reported_ties)",
        bound = "lower"
      ),
      m_upper = marginal_effect_bound(
        model = i,
        treatment_name = "treat_rank",
        moderator_name = "factor(rev_rank_reported_ties)",
        bound = "upper"
      )
    ) |> 
      filter(m != "2")) # remove Peru-Russia (not ranked 1,2,3)

mechanism_plots_df <- bind_rows(mechanism_plots_df, .id = "Outcome") |> 
  mutate(Outcome = factor(Outcome,
                          levels = c("Stop COVID-19 spread",
                                     "Help respondent country",
                                     "Increase support for sender",
                                     "Increase dependence on sender",
                                     "Obtain economic profits"
                          )))

mechanism_plots <- lapply(unique(mechanism_plots_df$Outcome), function(i)
  
  mechanism_plots_df |> 
    filter(Outcome == i) |> 
  ggplot(aes(x = m,
             y = m_effect)) +
    geom_point() +
    geom_errorbar(aes(x = m, ymin = m_lower, ymax = m_upper, width = 0)) +
    geom_hline(yintercept = 0, linetype = "dotted") +
    scale_x_discrete("Reversed rank",
                     labels = c("3", "4", "5")) +
    scale_y_continuous("Marginal effect of treatment") +
    ggtitle(i) +
    theme(axis.line = element_line(size = 0.5),
          axis.text = element_text(size = 17),
          axis.title = element_text(size = 20),
          plot.title = element_text(size = 20, hjust = 0.5),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank())
  
)

names(mechanism_plots) <- c("stopcv",
                            "help",
                            "supp",
                            "dep",
                            "econ")

for (i in names(mechanism_plots)) {
  
  ggsave(filename = paste0("Figures/mechanism_marginaleffects_", i, ".pdf"),
         plot = mechanism_plots[[i]],
         width = 1*5,
         height = 1*5,
         units = "in",
         dpi = 320)
}

rm(diplomacy_long,
   mechanism_models,
   mechanism_plots,
   mechanism_plots_df,
   rank_model,
   rank_model_binary,
   rank_plot_density,
   rank_plots,
   rank_plots_df,
   i,
   mechanism_outcomes,
   marginal_effect,
   marginal_effect_bound,
   marginal_effects_models)

# Figure A1: Distribution of trust in vaccine-developing countries ----

diplomacy_wide <- read.csv("vaccine_diplomacy_wide.csv",
                           encoding = "UTF-8")

diplomacy_wide <- diplomacy_wide |>
  mutate(
    trust_china_rec = case_when(
      trust_china_end_pre==1 ~ "No Trust", 
      trust_china_end_pre==2 ~ "Little Trust", 
      trust_china_end_pre==3 ~ "Some Trust", 
      trust_china_end_pre==4 ~ "A lot of trust", 
      trust_china_end_pre==2.5 ~  "Don't Know"
    ),
    trust_india_rec = case_when(
      trust_india_end_pre==1 ~ "No Trust", 
      trust_india_end_pre==2 ~ "Little Trust", 
      trust_india_end_pre==3 ~ "Some Trust", 
      trust_india_end_pre==4 ~ "A lot of trust", 
      trust_india_end_pre==2.5 ~  "Don't Know"
    ),
    trust_russia_rec = case_when(
      trust_russia_end_pre==1 ~ "No Trust", 
      trust_russia_end_pre==2 ~ "Little Trust", 
      trust_russia_end_pre==3 ~ "Some Trust", 
      trust_russia_end_pre==4 ~ "A lot of trust", 
      trust_russia_end_pre==2.5 ~  "Don't Know"
    ),
    trust_uk_rec = case_when(
      trust_uk_end_pre == 1 ~ "No Trust",
      trust_uk_end_pre == 2 ~ "Little Trust",
      trust_uk_end_pre == 3 ~ "Some Trust",
      trust_uk_end_pre == 4 ~ "A lot of trust",
      trust_uk_end_pre == 2.5 ~  "Don't Know"
    ),
    trust_us_rec = case_when(
      trust_us_end_pre == 1 ~ "No Trust",
      trust_us_end_pre == 2 ~ "Little Trust",
      trust_us_end_pre == 3 ~ "Some Trust",
      trust_us_end_pre == 4 ~ "A lot of trust",
      trust_us_end_pre == 2.5 ~ "Don't Know"
    )
  ) |> 
  mutate(across(
    c("trust_china_rec",
      "trust_india_rec",
      "trust_russia_rec",
      "trust_uk_rec",
      "trust_us_rec"),
    ~ factor(., levels = c("No Trust", "Little Trust", "Don't Know", "Some Trust", "A lot of trust"))))

trust_df <- data.frame(trust = factor(c("No Trust", "Little Trust", "Don't Know", "Some Trust", "A lot of trust"), levels = c("No Trust", "Little Trust", "Don't Know", "Some Trust", "A lot of trust")))

trust_df <- trust_df |> 
  mutate(trust_china_prop = prop.table(table(diplomacy_wide$trust_china_rec)),
         trust_india_prop = prop.table(table(diplomacy_wide$trust_india_rec)),
         trust_russia_prop = prop.table(table(diplomacy_wide$trust_russia_rec)),
         trust_uk_prop = prop.table(table(diplomacy_wide$trust_uk_rec)),
         trust_us_prop = prop.table(table(diplomacy_wide$trust_us_rec)))

trust_df <- pivot_longer(trust_df, 
                     cols = trust_china_prop:trust_us_prop,
                     names_to = "country",
                     values_to = "prop") |> 
  mutate(country = case_when(
    country == "trust_china_prop" ~ "China",
    country == "trust_india_prop" ~ "India",
    country == "trust_russia_prop" ~ "Russia",
    country == "trust_uk_prop" ~ "United Kingdom",
    country == "trust_us_prop" ~ "United States"),
    prop = as.numeric(prop))

trust_plots <- lapply(sort(unique(trust_df$country)), function(i)
  
  trust_df |> 
    filter(country == i) |> 
    ggplot(aes(x=factor(trust), y = prop, fill=factor(trust))) +
    geom_bar(stat = 'identity', color = 'black', fill='lightblue') +
    ylab("") +
    xlab("") +
    ggtitle(i) +
    scale_y_continuous(labels = scales::percent_format(accuracy = 5L)) +
    scale_x_discrete(labels = c("No\n Trust", 
                                "Little\n Trust", 
                                "Don't\n Know", 
                                "Some\n Trust", 
                                "A lot of\n Trust")) +
    theme(plot.title = element_text(size = 23, hjust = 0.5)) +
    theme(axis.text = element_text(size=17)) +
    coord_cartesian(ylim=c(0, .4)) 
                      
)

names(trust_plots) <- paste0(c("china", "india", "russia", "uk", "us"), "_trust_outcome.pdf")

trust_plots

for (i in names(trust_plots)) {
    
    ggsave(filename = paste0("Figures/", i),
           plot = trust_plots[[i]],
           width = 1*5,
           height = 1*5,
           units = "in",
           dpi = 320)
  
}

rm(diplomacy_wide, trust_df, trust_plots, i)